# load in libraries
from equadratures import *
import numpy as np
import matplotlib.pyplot as plt
import plotly.graph_objects as go

pressure assumption: $\frac{p}{p_0} = y = R^{-3}f_1$ ^[1]\ density assumption: $\frac{\rho}{\rho_0} = \psi$ ^[2]\ velocity assumption: $u = R^{-3/2}\phi_1$ ^[3]\ where $f_1$, $\psi$, and $\phi_1$ are functions of $\eta = \frac{r}{R}$ where $r$ is the radial position and $R$ is the radius of the shock wave

Using the Navier-Stokes equation: $\rho \frac{Du}{Dt} = -\nabla p + \rho g + \mu \nabla^2 u$ under the assumption of no gravitational or viscous effects we have $\rho \frac{Du}{Dt} = -\nabla p$, using the defintion of a materal derivative and the equations for denisty, pressure, and velocity from above we get ^[4]:

$$ \large \rho \left(\frac{\partial u}{\partial t} + (u \cdot \nabla)u\right) = -\frac{\partial p}{\partial r} \longrightarrow \rho \left(\frac{\partial u}{\partial t} + (u \cdot \nabla)u\right) = -\frac{p_0}{\rho}\frac{\partial y}{\partial r} $$$$ \large \left(\frac{\partial (R^{-3/2}\phi_1)}{\partial t} + R^{-3/2}\phi_1 \left(\frac{\partial (R^{-3/2} \phi_1)}{\partial r}\right)\right) = -\frac{p_0}{\rho_0\psi} \frac{\partial ( R^{-3} f_1)}{\partial r} $$$$ \large R^{-3/2}\frac{\partial \phi_1}{\partial t} - \frac{3}{2}R^{-5/2}\frac{dR}{dt}\phi_1 + R^{-3/2}\phi_1 \left(R^{-3/2}\frac{\partial \phi_1}{\partial r}\right) = -\frac{p_0}{\rho_0 \psi} R^{-3} \frac{\partial f_1}{\partial r} $$

Since $f_1$ and $\phi_1$ are functions of $\eta$ their partial derivatives can be expanded:

$$ \large R^{-3/2}\frac{\partial \phi_1}{\partial \eta}\frac{\partial \eta}{\partial t}- \frac{3}{2}R^{-5/2}\frac{dR}{dt}\phi_1 + R^{-3/2}\phi_1 \left(R^{-3/2}\frac{\partial \phi_1}{\partial \eta}\frac{\partial \eta}{\partial r}\right) = -\frac{p_0}{\rho_0 \psi} R^{-3} \frac{\partial f_1}{\partial \eta}\frac{\partial \eta}{\partial r} $$

Writing $\frac{\partial \phi_1}{\partial \eta}$ as $\phi_1'$ and $\frac{\partial f_1}{\partial \eta}$ as $f_1'$ and using $\frac{\partial \eta}{\partial t} = \frac{\partial r/R}{\partial t} = -\frac{r}{R^2} \frac{dR}{dt}$ and $\frac{\partial \eta}{\partial r}=\frac{\partial r/R}{\partial t}=\frac{1}{R}$ we get that:

$$ \large -R^{-3/2}\phi_1'\frac{r}{R^2}\frac{dR}{dt}- \frac{3}{2}R^{-5/2}\frac{dR}{dt}\phi_1 + R^{-3/2}\phi_1 \left(R^{-3/2}\phi_1'\frac{1}{R}\right) = -\frac{p_0}{\rho_0 \psi} R^{-3} f_1'\frac{1}{R} $$

Simplifying this we have ^[5]:

$$ \large -\left(\frac{3}{2}\phi_1 + \eta\phi_1'\right)R^{-5/2}\frac{dR}{dt} + R^{-4}\left(\phi_1\phi_1' + \frac{p_0}{\rho_0 \psi}f_1'\right) = 0 $$

This equation can be satisfied when $\frac{dR}{dt}=AR^{-3/2}$ ^[6] where A is some constant and in the end we could simplify ^[5] to get ^[7]):

$$ \large -\left(\frac{3}{2}\phi_1 + \eta\phi_1'\right)R^{-5/2}AR^{-3/2} + R^{-4}\left(\phi_1\phi_1' + \frac{p_0}{\rho_0 \psi}f_1'\right) = 0 $$$$ \large \longrightarrow -A\left(\frac{3}{2}\phi_1 + \eta\phi_1'\right) + \left(\phi_1\phi_1' + \frac{p_0}{\rho_0 \psi}f_1'\right) = 0 $$

Continuing on and using the equation of continuity $\frac{\partial \rho}{\partial t} + u \frac{\partial \rho}{\partial r} + \rho\left(\frac{\partial u}{\partial r} + \frac{2u}{r}\right)=0$ ^[8] we can use the same equations for denisty and velocity:

$$ \large \rho_0 \frac{\partial \psi}{\partial t} + \rho_0R^{-3/2}\phi_1 \frac{\partial \psi}{\partial r} + \rho_0\psi \left(R^{-3/2}\frac{\partial \phi_1}{\partial r} + \frac{2R^{-3/2}\phi_1}{r}\right)=0 $$$$ \large \rho_0 \frac{\partial \psi}{\partial \eta}\frac{\partial \eta}{\partial t} + \rho_0R^{-3/2}\phi_1 \frac{\partial \psi}{\partial \eta}\frac{\partial \eta}{\partial r} + \rho_0\psi \left(R^{-3/2}\frac{\partial \phi_1}{\partial \eta}{\partial \eta}{\partial r} + \frac{2R^{-3/2}\phi_1}{r}\right)=0 $$$$ \large \rho_0 \psi'\left(-\frac{r}{R^2}\frac{dR}{dt}\right)+ \rho_0R^{-3/2}\phi_1 \psi'\frac{1}{R} + \rho_0\psi \left(R^{-3/2}\phi_1'\frac{1}{R} + \frac{2R^{-3/2}\phi_1}{r}\right)=0 $$$$ \large \rho_0 \psi'\left(-\frac{r}{R^2}AR^{-3/2}\right)+ \rho_0R^{-5/2}\phi_1 \psi' + \rho_0\psi \left(R^{-5/2}\phi_1' + \frac{2R^{-3/2}\phi_1}{r}\right)=0 $$$$ \large -A\eta \psi'R^{-5/2}+ R^{-5/2}\phi_1 \psi' + \psi \left(R^{-5/2}\phi_1' + \frac{2R^{-5/2}\phi_1}{\eta}\right)=0 $$

Simplifying this we have ^[9]:

$$ \large -A\eta \psi'+ \phi_1 \psi' + \psi \left(\phi_1' + \frac{2\phi_1}{\eta}\right)=0 \ (9) $$

For a perfect gas $p\rho^{-\gamma} = const$ therefore $\frac{D(p\rho^{-\gamma})}{Dt}=0$ and therefore $\frac{\partial (p\rho^{-\gamma})}{\partial t} + u \frac{\partial (p \rho^{-\gamma})}{\partial r} =0$ ^[10]

$$ \large p \frac{\partial \rho^{-\gamma}}{\partial t} + \rho^{-\gamma} \frac{\partial p}{\partial t} + u\left(p \frac{\partial \rho^{-\gamma}}{\partial r} + \rho^{-\gamma} \frac{\partial p}{\partial r}\right)=0 $$$$ \large -\gamma p \rho^{-\gamma-1}\frac{\partial \rho}{\partial t} + \rho^{-\gamma} \frac{\partial p}{\partial t} + u\left(-\gamma p \rho^{-\gamma-1}\frac{\partial \rho}{\partial r} + \rho^{-\gamma} \frac{\partial p}{\partial r}\right)=0 $$$$ \large -\gamma p \rho^{-1}\frac{\partial \rho}{\partial t} + \frac{\partial p}{\partial t} + u\left(-\gamma p \rho^{-1}\frac{\partial \rho}{\partial r} + \frac{\partial p}{\partial r}\right)=0 $$$$ \large -\gamma p_0R^{-3}f_1 \rho_0(\rho_0 \psi)^{-1}\frac{\partial \psi}{\partial t} -3 f_1 p_0 R^{-4}\frac{dR}{dt} + p_0R^{-3}\frac{\partial f_1}{\partial t} + u\left(-\gamma p_0R^{-3}f_1 \rho_0(\rho_0 \psi)^{-1}\frac{\partial \psi}{\partial r} + p_0R^{-3}\frac{\partial f_1}{\partial r}\right)=0 $$$$ \large -\gamma p_0R^{-3}f_1 \psi^{-1} \psi' \left(-\frac{r}{R^2}\frac{dR}{dt}\right) -3 f_1 p_0 R^{-4}\frac{dR}{dt} + R^{-3}p_0f_1' \left(-\frac{r}{R^2}\frac{dR}{dt}\right) + u\left(-\gamma p_0R^{-3}f_1 \psi^{-1}\psi'\frac{1}{R} + p_0R^{-3}f_1'\frac{1}{R}\right)=0 $$$$ \large A\gamma p_0R^{-11/2}f_1\eta \psi^{-1} \psi' -3 f_1 p_0 R^{-4}AR^{-3/2} + A\eta R^{-11/2}p_0f_1' + R^{-3/2}\phi_1 \left(-\gamma p_0R^{-4}f_1 \psi^{-1}\psi' + p_0R^{-4}f_1'\right)=0 $$$$ \large A\gamma f_1\eta \psi^{-1} \psi' -3f_1A + A\eta f_1' + \phi_1 \left(-\gamma f_1 \psi^{-1}\psi' + f_1'\right)=0 $$$$ \large -A(3f_1 + \eta f_1') +\gamma f_1\psi^{-1} \psi'( A\eta -\phi_1) + \phi_1f_1'=0 $$

Simplifying this becomes ^[11]: $$ \large A(3f_1 + \eta f_1') +\gamma f_1\frac{\psi'}{\psi}(\phi_1 - A\eta) - \phi_1f_1'=0 $$

We can reduce ^[7], ^[9], and ^[11] by forming nondimensional forms of $f_1$ and $\phi_1$ where $f=f_1\frac{a^2}{A^2}$ (12) and $\phi = \frac{\phi_1}{A}$ (13) where $a^2 = \gamma \frac{p_0}{\rho_0}$ and is the speed of sound. Doing this for ^[7] we get ^[7a]:

$$ \large -A\left(\frac{3}{2}A\phi + \eta A\phi'\right) + \left(A^2\phi\phi' + \frac{p_0}{\rho_0 \psi}\frac{A^2}{a^2}f'\right) = 0 $$$$ \large -\frac{3}{2}\phi + \frac{p_0}{\rho_0 \psi}\frac{\rho_0}{\gamma p_0}f' = \eta \phi' - \phi\phi' $$$$ \large \phi'(\eta - \phi)= \frac{1}{\gamma}\frac{f'}{\psi} -\frac{3}{2}\phi $$

Using the subsitutions in ^[9] gives ^[9a]:

$$ \large -A\eta \psi'+ A\phi \psi' + \psi \left(A\phi' + \frac{2A\phi}{\eta}\right)=0 $$$$ \large \psi \left(\phi' + \frac{2\phi}{\eta}\right)=\psi'(\eta -\phi) $$$$ \large \frac{\psi'}{\psi} = \frac{\left(\phi' + \frac{2\phi}{\eta}\right)}{\eta -\phi} $$

Using the substitution in ^[11] gives ^[11a]:

$$ \large A(3\frac{A^2}{a^2}f + \eta \frac{A^2}{a^2}f') +\gamma \frac{A^2}{a^2}f\frac{\psi'}{\psi}(A\phi - A\eta) - A\phi\frac{A^2}{a^2}f'=0 $$$$ \large 3f + \eta f' +\gamma f\frac{\psi'}{\psi}(\phi - \eta) - \phi f'=0 $$

Substituting ^[7a] and ^[9a] into ^[11a] gives ^[14]:

$$ \large 3f + \eta f' +\gamma f\frac{\left(\phi' + \frac{2\phi}{\eta}\right)}{\eta -\phi}(\phi - \eta) - \phi f'=0 $$$$ \large 3f + \eta f' -\gamma f\left(\phi' + \frac{2\phi}{\eta}\right) - \phi f'=0 $$$$ \large 3f + \eta f' -\gamma f\left(\frac{\frac{1}{\gamma}\frac{f'}{\psi} - \frac{3}{2}\phi}{\eta - \phi} + \frac{2\phi}{\eta}\right) - \phi f'=0 $$$$ \large f'\left((\eta - \phi)^2 - \frac{f}{\psi}\right) = f\left(-3\eta + \phi\left(3 + \frac{1}{2}\gamma\right) - \frac{2\gamma \phi^2}{\eta}\right) $$

Using the Rankine-Hugoniot relations and reducing we get ^[15], ^[16], and ^[17]:

$$ \large \frac{\rho_1}{\rho_0} = \frac{\gamma - 1 +(\gamma+1)y_1}{\gamma+1 + (\gamma-1)y_1} $$$$ \large \frac{U^2}{a^2} = \frac{1}{2\gamma}(\gamma -1 + (\gamma+1)y_1) $$$$ \large \frac{u_1}{U} = \frac{2(y_1-1)}{\gamma-1 + (\gamma+1)y_1} $$

Where $\rho_1$, $y_1$, and $u_1$ are the values of $\rho$, $\gamma$, and $y_1$ immediately behind/before the shock and $U=\frac{dR}{dt}$, we can simplify these for the case where $y_1 >>$ to get ^[15a], ^[16a], and ^[17a]:

$$ \large \frac{\rho_1}{\rho_0} = \frac{\gamma+1}{\gamma-1} $$$$ \large \frac{U^2}{a^2} = \frac{\gamma +1}{2\gamma}y_1 $$$$ \large \frac{u_1}{U} = \frac{2}{\gamma+1} $$

Since $\psi=\frac{\rho}{\rho_0}$ ^[15a] becomes ^[15b]:

$$ \large \psi = \frac{\gamma+1}{\gamma-1} $$

Next, since $f=\frac{a^2}{A^2}$ then $a^2=\frac{A^2f}{f_1}$ and because $U=\frac{dR}{dt}=AR^{-3/2}$ then $U^2=A^2R^{-3}$, additionally $y=R^{-3}f_1$ using these equations ^[16a] becomes ^[16b]:

$$ \large \frac{A^2R^{-3}f_1}{A^2f} = \frac{\gamma + 1}{2\gamma}R^{-3}f_1 $$$$ \large f = \frac{2\gamma}{\gamma + 1} $$

Lastly, for ^[17a] we can use that $u=R^{-3/2}A\phi$ and $U=\frac{dR}{dt}=AR^{-3/2}$ to get ^[17b]:

$$ \large \frac{AR^{-3/2}\phi}{AR^{-3/2}} = \frac{2}{\gamma+1} $$$$ \large \phi = \frac{2}{\gamma+1} $$

Note: In Taylor's paper ^[16a] is different: However, this does not derive from equation 16 and does not simplify to ^[16b] so I'm not sure what he did but the end result of ^[16b] is the same

By using (^[15b], ^[16b], ^[17b]) and using that $\gamma=1.4$ we can get values for $f$, $\phi$, and $\psi$ at $\eta=1$:

# define values at eta=1
f_0 = 7/6
phi_0 = 5/6
psi_0 = 6
g = 1.4

Using ^[14], ^[7a], and ^[9a] we can define functions to calculation $f'$, $\phi'$, and $\psi'$

# define f'
def dev_f(f, phi, psi, eta):
    f_prime = (f*(-3*eta+phi*(3+0.5*g)-(2*g*phi**2)/eta))/((eta-phi)**2 -f/psi)
    return f_prime 

# define phi'
def dev_phi(f, phi, psi, eta):
    f_prime = dev_f(f, phi, psi, eta)
    phi_prime = ((1/g)*(f_prime/psi)-(3/2)*(phi))/(eta-phi)
    return phi_prime

# define psi'
def dev_psi(f, phi, psi, eta):
    phi_prime = dev_phi(f, phi, psi, eta)
    psi_prime = psi*(phi_prime+2*phi/eta)/(eta-phi)
    return psi_prime

To find an approximate function for $f$, $\phi$, and $\psi$ we can calcuate the change at very small increments

etas = np.linspace(1,0.5, 100)

# step by step calculation 
fs = []
phis = []
psis = []

f = f_0
phi = phi_0
psi = psi_0

f_prime_list = []

delta = np.mean(np.diff(etas))
for eta in etas:
    fs.append(f)
    phis.append(phi)
    psis.append(psi)
    
    f_prime = dev_f(f, phi, psi, eta)
    phi_prime = dev_phi(f, phi, psi, eta)
    psi_prime = dev_psi(f, phi, psi, eta) 
    f_prime_list.append(f_prime)
    
    f += delta*f_prime
    phi += delta*phi_prime
    psi += delta*psi_prime

The output for $\phi$ is very nearly linear and can be approximated using ^[14] and ^[9a]:

$$ \large 3f + \eta f' -\gamma f\left(\phi' + \frac{2\phi}{\eta}\right) - \phi f'=0 $$

Note: in Taylor ^[19] is written differently but both forms are equivalent \ This can be simplifed to get ^[20]

$$ \large \frac{f'}{f}(\eta - \phi) = \gamma \phi' + \frac{2\gamma \phi}{\eta} - 3 $$

If the left-hand side is ignored then the solution is $\phi = \frac{\eta}{\gamma}$ ^[21] this linear relationship follows what is calculated by the step process however it starts to break away at high values of $\eta$ meaning there is a correction factor of the form $\alpha\eta^n$ where $\alpha$ and n are constants ^[22]

$$ \large \phi = \frac{\eta}{\gamma} + \alpha\eta^n $$

These constants can be found using ^[16b] and ^[17b], firstly by using ^[17b] and plugging that value in for $\phi$ at $\eta=1$ we get ^[23]:

$$ \large \frac{2}{\gamma +1} = \frac{1}{\gamma} + \alpha (1)^n \longrightarrow \alpha = \frac{\gamma -1}{\gamma(\gamma +1)} $$

Next, we can use (^[15b], ^[16b], ^[17b]) to find the true value of $f'$ at $\eta=1$ getting that $f' = \frac{2\gamma^2 +7\gamma-3}{\gamma-1}$, next we can take the derivative of ^[22] to get $\phi' = \frac{1}{\gamma} +\alpha n \eta^{n-1}$ then plug this value into ^[20] and set them equal in order to calculate n ^[24]:

$$ \large n = \frac{7\gamma-1}{\gamma^2 -1} $$

Then we can use ^[22] to solve for $f$ and $\psi$, firstly plugging ^[22] into ^[20] we get ^[25]:

$$ \large \frac{f'}{f} = \frac{(n+2)\alpha\gamma^2\eta^{n-2}}{\gamma-1-\gamma\alpha\eta^{n-1}} $$

Taking the derivative and solving for the derivation constant using the value of $f$ at $\eta=1$ we get that ^[26]:

$$ \large ln(f) = ln\left(\frac{2\gamma}{\gamma+1}\right) - \frac{2\gamma^2+7\gamma-3}{7-\gamma}ln\left(\frac{\gamma+1-\eta^{n-1})}{\gamma}\right) $$

Next we can use ^[22] and its derivative to solve for $\psi$ by using ^[9a]:

$$ \large \frac{\psi'}{\psi} = \frac{\frac{1}{\gamma} + \alpha n \eta^{n-1} + \frac{2}{\eta}\left(\frac{\eta}{\gamma} + \alpha \eta^n\right)}{\eta - \left(\frac{\eta}{\gamma} + \alpha\eta^n\right)}=\frac{3+(n+2)\gamma\alpha\eta^{n-1}}{(\gamma-1)\eta-\alpha\gamma\eta^n} $$

Taking the derivative and solving for $\psi$ gets ^[28]:

$$ \large ln(\psi) = ln\left(\frac{\gamma+1}{\gamma-1}\right)+\frac{3}{\gamma-1}ln(\eta)-2\frac{\gamma+5}{7-\gamma}ln\left(\frac{\gamma+1-\eta^{n-1}}{\gamma}\right) $$

Below we can check these equations for $\phi$, $f$, and $\psi$

# guess constants and values
a = (g-1)/(g*(g+1))
n = (7*g-1)/(g**2-1)
phi_guess = etas/g
phi_guess_corrected = etas/g + a*etas**n
f_guess_corrected = (2*g/(g+1))*(((g+1)/g)-((etas**(n-1))/g))**(-(2*g**2+7*g-3)/(7-g))
psi_guess_corrected = ((g+1)/(g-1))*(etas**(3/(g-1)))*((g+1-etas**(n-1))/g)**(-2*(g+5)/(7-g))

# taylors data
taylor_etas = np.linspace(1,0.5,26)
taylor_fs = [f_0, 0.949, 0.808, 0.711, 0.643, 0.593, 0.556, 0.528, 0.507, 0.491, 0.478, 0.468, 0.461, 0.455, 0.45, 0.447, 
             0.444, 0.442, 0.44, 0.439, 0.438, 0.438, 0.437, 0.437, 0.437, 0.436]
taylor_phis = [phi_0, 0.798, 0.767, 0.737, 0.711, 0.687, 0.665, 0.644, 0.625, 0.607, 0.590, 0.573, 0.557, 0.542, 0.527, 
               0.513, 0.498, 0.484, 0.47, 0.456, 0.443, 0.428, 0.415, 0.402, 0.389, 0.375]
taylor_psis = [psi_0, 4, 2.808, 2.052, 1.534, 1.177, 0.919, 0.727, 0.578, 0.462, 0.37, 0.297, 0.239, 0.191, 0.152, 0.12, 
               0.095, 0.074, 0.058, 0.044, 0.034, 0.026, 0.019, 0.014, 0.01, 0.007]

fig = go.Figure()

fig.add_scatter(x=etas, y=fs, name='step', line_color = 'skyblue')
fig.add_scatter(x=taylor_etas, y=taylor_fs, name='Taylor', line_color='darkviolet', line_dash='dash')
fig.add_scatter(x=etas, y=f_guess_corrected, name='corrected guess', line_dash='dot', line_color='deeppink')

fig.add_scatter(x=etas, y=phis, visible=False, name='step', line_color='skyblue')
fig.add_scatter(x=taylor_etas, y=taylor_phis, visible=False, name='Taylor', line_dash='dash', line_color='darkviolet')
fig.add_scatter(x=etas, y=phi_guess, visible=False, name='guess', line_dash='dashdot', line_color='pink')
fig.add_scatter(x=etas, y=phi_guess_corrected, visible=False, name='corrected guess', line_dash='dot', line_color='deeppink')

fig.add_scatter(x=etas, y=psis, visible=False, name='step', line_color='skyblue')
fig.add_scatter(x=taylor_etas, y=taylor_psis, visible=False, name='Taylor', line_dash='dash', line_color='darkviolet')
fig.add_scatter(x=etas, y=psi_guess_corrected, visible=False, name='corrected guess', line_dash='dot', line_color='deeppink')

fig.update_layout(updatemenus=[dict(active=0, buttons=[dict(label='f', method='update', 
                  args=[{'visible':[True, True, True, False, False, False, False, False, False, False]}]), 
                  dict(label=r'$\phi$', method='update', args=[{'visible':[False, False, False, True, True, True, True, 
                  False, False, False]}]), dict(label=r'$\psi$', method='update', args=[{'visible':[False, False, False, 
                  False, False, False, False, True, True, True]}])])])

fig.update_xaxes(title=r'$\eta \left(\frac{r}{R}\right)$')

fig.show()

To calculate the energy of the explosion we need to look at two parts, kinetic and heat \ Kinetic energy in integral form using spherical coordinates is $KE=4\pi\int_0^Rfrac{1}{2}\rho u^2r^2dr$ \ In terms of the variables $\phi$, $\psi$, $f$, and $\eta$ we can write this as:

$$ \large 2\pi\int_0^11 \rho_0\psi A^2R^{-3}\phi^2 R^2\eta^2 Rd\eta = 2\pi A^2\rho_0\int_0^1 \psi \phi^2 \eta^2 d\eta $$

Heat energy in integral form using spherical coordinates is $HE=4\pi\int_0^R\frac{pr^2}{\gamma-1}dr$ \ In terms of the variables $\phi$, $\psi$, $f$, and $\eta$ we can write this as:

$$ \large \frac{4\pi}{\gamma-1} \int_0^1 p_0 \frac{A^2}{a^2}fR^{-3} R^2\eta^2 Rd\eta = \frac{4\pi A^2 \rho_0}{\gamma(\gamma-1)} \int_0^1 f \eta^2 d\eta $$

Therefore the total Energy of the system can be written as:

$$ \large E = 2\pi A^2 \rho_0 \left(\int_0^1 \psi \phi^2 \eta^2 d\eta + \frac{2}{\gamma(\gamma-1)}\int_0^1 f \eta^2 d\eta\right) $$

Using that $\frac{dR}{dt}=AR^{-3/2}$ we can integrate and solve for A to get that $A=\frac{2}{5}R^{5/2}t^{-1}$ therefore:

$$ \large E = 2\pi \frac{4}{25}R^5t^{-2} \rho_0 \left(\int_0^1 \psi \phi^2 \eta^2 d\eta + \frac{2}{\gamma(\gamma-1)}\int_0^1 f \eta^2 d\eta\right) $$

We can state that $K=\frac{4}{25}\left(2\pi\int_0^1 \psi \phi^2 \eta^2 d\eta + \frac{4\pi}{\gamma(\gamma-1)}\int_0^1 f \eta^2 d\eta\right)$ so that the energy is $E=K\rho_0R^5t^{-2}$

eta_p = Parameter(lower=0, upper=1, order=10)
basis = Basis('tensor-grid')
poly = Poly(parameters=eta_p, basis=basis, method='numerical-integration')
pts, wts = poly.get_points_and_weights()

def K(gamma, eta):
    a = (gamma-1)/(gamma*(gamma+1))
    n = (7*gamma - 1)/(gamma**2 - 1)
    phi = (eta/gamma) + a*eta**n
    log_psi = np.log((gamma+1)/(gamma-1)) + 3/(gamma-1) * np.log(eta) - 2*(gamma+5)/(7-gamma) * \
              np.log((gamma+1-eta**(n-1))/gamma)
    psi = np.exp(log_psi)
    log_f = np.log(2*gamma/(gamma+1)) - (2*gamma**2 + 7*gamma - 3)/(7-gamma) * np.log((gamma+1-eta**(n-1))/gamma)
    f = np.exp(log_f)
    I_1 = wts @ (eta**2 * phi**2 * psi)
    I_2 = wts @ (eta**2 * f)
    return 4/25 * (2*np.pi*I_1 + 4*np.pi/(gamma*(gamma-1)) * I_2)

gammas = np.linspace(1.2, 1.7, 100)

Ks = []
for gamma in gammas:
    Ks.append(K(gamma, pts))

rho_0 = 1.25 # kg/m^3
ergs = 1e7 # ergs/J
tnt = 4.25e16 # ergs/tons TNT

From that $\frac{dR}{dt}=AR^{-3/2}$ the relationship between R and t can be found by taking the integral: $At=\frac{2}{5}R^{5/2}$ meaning that the time is proportional to some constant times $R^{5/2}$ with the values taken from the images of the explosion we can find this constant:

R = np.array([11.1, 19.9, 25.4, 28.8, 31.9, 34.2, 36.3, 38.9, 41, 42.8, 44.4, 46, 46.9, 48.7, 59, 61.1, 62.9, 64.3, 65.6,
             67.3, 106.5, 130, 145, 175, 185]) #m
t = np.array([0.1, 0.24, 0.38, 0.52, 0.66, 0.8, 0.94, 1.06, 1.22, 1.36, 1.5, 1.65, 1.79, 1.93, 3.26, 3.53, 3.8, 4.07, 4.34,
              4.61, 15, 25, 34, 53, 62])*(1e-3) #millisec -> sec

# find the approximate value of R^5 t^-2
R5t2 = np.mean((R**5)*(t**(-2)))

# find the total energy based on gamma and convert to tons TNT
E = np.array(Ks)*rho_0*R5t2
E_ergs = E*ergs
E_TNT = E_ergs / tnt
E_taylor = [14.4e20, 9.74e20, 7.14e20, 4.06e20] # ergs

# plot the energy vs the values that Taylor got 
fig = go.Figure()

fig.add_scatter(x=gammas, y=E_TNT[:,0], line_color='black')
fig.add_hline(y=20000, line_dash='dash', annotation_text='Actual Estimated Energy of Explosion')
fig.add_scatter(x=[1.2, 1.3, 1.4, 1.667], y=E_taylor, yaxis='y2', mode='markers', marker_color='blue',
                name='Taylor')
fig.add_scatter(x=gammas, y=E_ergs[:,0], yaxis='y2', line_color='black', name='Code')

fig.update_layout(yaxis=dict(title='Tons TNT', side='right'), yaxis2=dict(title='Energy (ergs)'), showlegend=False,
                  xaxis=dict(range=[1.15,1.75]))

fig.show()

First test released an approximate amount of 20,000 tons of TNT (source: Wikipedia. ) )

The pressure ratio is $\frac{p}{p_0}=R^{-3}f_1=\frac{A^2}{a^2}R^{-3}f$ using that $a^2=\gamma\frac{p_0}{\rho_1}$ we can get that $p=p_0A^2\frac{\rho_0}{\gamma p_0}R^{-3}f$, next using that $A=\frac{2}{5}R^{5/2}t^{-1}$ the pressure is $p=\frac{\rho_0}{\gamma}\frac{4}{25}R^{5}t^{-2}R^{-3}f$ to find the max pressure at each value of gamma we need the highest value of f which happens at $\eta=1$ and is $f=\frac{2\gamma}{\gamma+1}$ so the pressure becomes $p_{max}=\frac{2\rho_0}{\gamma+1}\frac{4}{25}R^{5}t^{-2}R^{-3}$ \ The energy is equivalent to $E=K\rho_0 R^{5}t^{-2}$ so the pressure is $p_{max}=\frac{8}{25}\frac{E}{K(\gamma+1)}R^-3$ where $\frac{E}{K}$ is a constant value

p_z = []
for gamma in gammas:
    g_p = []
    for i, r in enumerate(R):
        p_max = (rho_0)*((r**5)*(t[i]**(-2)))*(8/25)*(1/(1+gamma))*r**(-3)
        g_p.append(p_max)
    p_z.append(g_p)

# plot the pressure vs. the radius of the shock and the gamma value
fig = go.Figure()

fig.add_contour(z=np.log10(p_z), x=R, y=gammas, contours_coloring='heatmap', line_width=0, colorscale='viridis',
                colorbar=dict(title='Max Pressure [log(P)]', titleside='right'))
fig.update_layout(xaxis=dict(title='R [m]'), yaxis=dict(title=r'$\gamma$'))

fig.show()

Energy and pressure evaluations are taken from Taylor Part I equations 34-38 and Taylor Part II equations 1-8 as well as tables 1-3